SPLITTING SCHEMES FOR HYPERBOLIC HEAT CONDUCTION 

EQUATION 

PETR N. VABISHCHEVICH* 

Abstract. Rapid processes of heat transfer are not described by the standard heat conduction 
equation. To take into account a finite velocity of heat transfer, we use the hyperbolic model of heat 
conduction, which is connected with the relaxation of heat fluxes. In this case, the mathematical 
model is based on a hyperbolic equation of second order or a system of equations for the temperature 
and heat fluxes. In this paper we construct for the hyperbolic heat conduction equation the additive 
schemes of splitting with respect to directions. Unconditional stability of locally one-dimensional 
splitting schemes is established. New splitting schemes are proposed and studied for a system of 
equations written in terms of the temperature and heat fluxes. 
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1. Introduction. Linear parabolic theories of diffusion and heat conduction 
are based on the Fick and Fourier laws, respectively, and predict an infinite speed of 
propagation [5lll2|. In this case, the amplitude of propagating perturbations decreases 
exponentially with the distance and the infinite speed of perturbations can often be 
ignored. Nevertheless, in many applied problems the wave nature of heat transfer 
should be taken into account. Since paper [7 , various corrections have been proposed 
for parabolic heat conduction models in order to eliminate the paradox of infinite 
speed of perturbation propagation [13j [TH [26] . 

The standard parabolic heat conduction model is based on the explicit repre- 
sentation of the heat flux through the temperature gradient. The hyperbolic heat 
conduction model includes an additional term with the time derivative for the heat 
flux which is proportional to the relaxation tensor. More general models (see |13j ) 
in addition includes the relaxation of the temperature gradient. Separate attention 
should be given to the hyperbolic model of convection-diffusion for moving media 

mm- 

Two possibilities can be highlighted in constructing computational methods for 
the approximate solution of hyperbolic heat transfer problems. The first is connected 
with the transition from a system of the first order evolutionary equations for the 
temperature and heat flux to a single hyperbolic equation of second order. In contrast 
to the standard parabolic equation of heat conduction there does present a term with 
the second time derivative. The second possibility is based on the usage of the initial 
scalar-vector system of equations. 

To solve approximately the boundary value problem for the hyperbolic heat con- 
duction equation, classical numerical methods can be used including finite-difference 
approximations in space, finite volume schemes or finite-element approximations. For 
instance, three-level difference schemes for the hyperbolic heat conduction equation 
are constructed in [4] [21] . To investigate the stability and convergence of difference 
schemes, the general theory of stability for operator-difference schemes is used in 
[T8l [T9] . Investigation of the stability on the basis of a priori estimates of the finite- 
difference solution for the model with the relaxation of the temperature gradients 
is given in [6j [35]. An analysis of possibilities to use the simplest schemes of first 
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and second order for integration in time is given in |16) for model one-dimensional 
problems of hyperbolic heat conduction. 

The system of equations governing thermal processes in terms of the tempera- 
ture and heat flux has a defined structure with conjugated to each other operators. 
Such a structure of the mathematical model makes possible to use this feature in 
the construction of computational algorithms [U 117) . In a number of papers (see, 
eg, [3 [27l [29l [34]) the hyperbolic nature of mathematical models with heat waves 
emphasizes in using traditional technologies of compressible media dynamics. Var- 
ious finite element methods are applied in papers [51 [5] for equations of hyperbolic 
convection-diffusion theory. 

Much attention is paid to the construction of the additive schemes (splitting 
schemes) for approximate solving initial-boundary value problems for multi-dimensional 
partial differential equations [33J . Transition to a sequence of more simple prob- 
lems allows to construct, for example, economical difference schemes - schemes based 
on the splitting with respect to spatial variables. In some cases it is reasonable to 
perform splitting with respect to subproblcms of different nature - splitting in physical 
processes. At present regionally-additive schemes (domain decomposition methods) 
are actively discussed [24] . These schemes are oriented to the construction of compu- 
tational algorithms for parallel computers. 

Additive difference schemes in general conditions of the splitting of the prob- 
lem operator into a sum of noncommutative non-selfadjoint operators are obtained 
in the most simple way for the case of two-component splitting. In this case for the 
evolutionary equation of first order the classical alternative direction schemes, factor- 
ized and predictor-corrector schemes are unconditionally stable at weak restrictions. 
A more complicated situation takes place in the case of multi-component splitting 
(splitting into three and more operators). For these problems the most interesting 
results are obtained on the basis of the concept of summarized approximation. The 
initial problem at the transition from one time level to another is divided into several 
subproblems, and each of these subproblems, in general, do not approximate the ini- 
tial problem. On this way, unconditionally stable schemes of componentwise splitting 
(locally one-dimensional schemes of splitting with respect to spatial variables) are 
constructed. 

A new class of operator-difference splitting schemes - vector additive schemes - 
was developed in papers [TJ [3D] . In this class of schemes we go from the initial scalar 
problem for one unknown function to the problem for a vector, each component of 
which can be treated as the solution of the problem. On this way we construct the 
full approximation schemes for evolutionary equations of the first and second order 
based on a general multi-component splitting. New additive difference schemes for 
differential-operator equations of the first and second order for the general case of 
splitting with an arbitrary number of pairwise noncommutative operator terms were 
constructed in |23[ 124] using the principle of regularization. 

The main theoretical results on the stability and convergence of the additive 
schemes were obtained for scalar evolutionary equations of the first order and, in 
some cases, for second-order equations. Splitting schemes for systems of evolution- 
ary equations are of considerable interest for computational practice. For standard 
parabolic and hyperbolic systems of equations with selfadjoint elliptic operators addi- 
tive schemes were constructed in |I8] using the principle of regularization for difference 
schemes. The Cauchy problem for a special linear system of first order equations in 
the Hilbert space with the conjugate operators (divergence and gradient) is considered 



SPLITTING SCHEMES FOR HYPERBOLIC HEAT 



3 



in paper |31j . Such a structure of equations is characteristic for the considering here 
problems of hyperbolic heat transfer. 

In the present work there are constructed splitting schemes with respect to spatial 
variables for the approximate solving the equation of hyperbolic heat conduction. 
Unconditionally stable locally one-dimensional difference schemes are constructed here 
both for a single heat conduction equation and for the system of equations based on 
the temperature and heat flux as unknowns. This paper is organized as follows. In 
section 2 the differential problem is formulated for the hyperbolic heat conduction. 
Appropriate a priori estimates are obtained for the solution of the hyperbolic equation 
in the both above mentioned formulations. Approximation in space is discussed in 
Section 3 for a model problem in a rectangle. It was shown that the grid operators of 
divergence and gradient are ajoint each other. Standard three-level difference schemes 
for the hyperbolic heat conduction equation are constructed in Section 4. The a priori 
estimates are derived for the difference solution. Difference schemes for the system 
of equations based on the temperature and heat flux as unknowns are considered 
in Section 5. Unconditionally stable schemes are derived via the regularization of 
explicit-implicit schemes. Locally one-dimensional schemes for the hyperbolic heat 
conduction equation are studied in Section 6. Splitting scheme for the system of 
hyperbolic heat conduction equations are proposed in Section 7. 

2. Differential problem. Temperature u(x, t) in bounded domain il with bound- 
ary dQ is governed by the equation 

du 

(2.1) c— +divq=/, xefJ, 0<<<T, 

ot 

where x = (x\,X2, ...,x n ) is a point in space, t is the time (T > 0), and q = q(x, t) 
is the heat flux. In ( ref (2.1)) c = c(x) > cq > is the specific heat capacity of 
a medium, and / = /(x, t) is the rate of volumetric heat sources. The standard 
(parabolic) model of the heat conduction results from the following representation for 
the heat flux (Fourier's law) 

(2.2) q + fcgradu = 0, 

where k = fc(x) > fco > is the thermal conductivity of the medium. Substitution of 
in (I2.1[) leads us to the parabolic heat conduction equation 



f)ii 

(2.3) c— - div(kgradu) = /, xe!l, < t < T, 

supplemented by appropriate boundary and initial conditions. 

In the model of the hyperbolic heat conduction instead of (|2.2I) we use the follow- 
ing relation 

do 

(2.4) q+z/-^ + fcgradu = 0, 

Ot 

where v is the relaxation parameter for the heat flux. From (|2.1j) and (|2.4j) we obtain 
the hyperbolic heat conduction equation 

Q 2 u du d f 

(2.5) vc— — +c— - div(fcgradw) =/ + v xett. < t < T. 
K ' dt 2 dt v 6 ' J dt 
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Consider a model boundary value problem for equation 
(EH])), where the boundary conditions are as follows 

(2.6) u(x, t)=0, x G dVl, 0<t<T. 

In addition, two initial conditions are prescribed 

du 



(system (|2.1I) and 



(2.7) 



u(x,0) = u (x), 



at 



(x,0) = ui(x), x e ft. 



The simplest a priori estimates for problem ([2"3 ]) -([27r ]) . ( (j2~T|) . (j2~4")l . ([2~p]) . ([277])) will 
be derived now in order to be our guidelines in the investigation of grid problems. 

Let (•, ■) and || • || be the scalar product and norm in "H = L2(fi), respectively. 
Multiplying scalar ly equation (I2.5[) by du/dt in H we obtain 



, du du 



2 dt 

The right hand side of 



du du 



dt' dt 



1 d 



2dt 



df du 



(2.9) 



df du\ 



is estimated as follows 



du du 



dt' dt J 

From (|2.8|) . ()2.9[) we have the inequality 

d 



4 IC 



5/ 



ft 



/ + "it ), f + V^I 



at' dt 



df 



dt 



(2.10) 
Here 
(2.11) 



,S < n 
dt ~ 2 



2\ c f + v W f + v m) 



, du du 



{k gradw, gradw) 



defines the squared norm for the solution of (|2.5[) - (I2.7[) with boundary conditions 
(|2.6[) . Applying to (|2 . 10[) the Gronwall lemma, we obtain the desired estimate 



(2.12) 



S(t) < S(0) + - 



,-1/2 



dO. 



At v — estimate (|2.12p degenerates into the corresponding estimate for the solution 
of parabolic heat equation (12.31) . 

For the system of equations instead of initial conditions ()2.7|) it is more natural 
to use 



(2.13) 



u(x, 0) = uo(x), q(x, 0) = go(x), x e dfl, 



ie instead of the rate of temperature variation we define the heat flux. The transition 
from one to another initial conditions is provided by equation (|2 . 3|) . 

To obtain a simple a priori estimate for system (|2.ip . (|2.4[) we scalar ly multiply 
equation (12. 1| by u, and (|2.4[) - by fc _1 q and sum them. This gives 



(2.14) 



1 d , \ v d ,. _i , ,, _■■ , 
(cm, u) + - — (fc q, q) + (fc q, q) 



2dt 



2 rfi 



(/,«). 
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For the right hand side we use the estimate 

(/, u )<i( cu , u ) + i( c -V,/). 

From (12.141) we obtain 

(2.15) | G < G+ ( c -i/, /); 

(2.16) G{t) = (cu,u) + iy(k- 1 q,q). 
From (|2.15p we derive the estimate 

t 

(2.17) G(t) < exp(t)G(Q) + /exp(i - 0)||c _1 / a /(x, 0)\\ 2 d6, 



which ensures the stability of the solution of system (|2.1j) . (|2.4p with respect to initial 
data (|2.13[) and the right hand side. 

3. Approximation in space. Let us consider the 2D model problem of the 
hyperbolic heat conduction in the rectangle 

Q = { x | x = (xi,X2), < x a < l a , a — 1,2}. 

Let q a , a = 1,2 be the Cartesian components of heat flux q = (91,52). The system 
of equations (|2.1[) . (|2.4I) in the coordinate- wise representation takes the form 

/q 1 \ du A dq a 

(3 - 1} c ^ + E^- = /> 

a— 1 

(3.2) q a + i/-Jp + k-^- = 0, a = 1,2. 

On the set of functions u, satisfying homogeneous boundary conditions (|2.6p . we 
define the operators 

(3.3) ^l Q u = — — , a = 1,2. 

ox a 

Taking into account that 

n n 

we have 

(3.4) A > = ~^ a = 1 ' 2 

for the conjugate operators. In view of (|3.3p . (|3.4I) the system of equations (|3.ip , (|3.2p 
with boundary conditions (|2.6p can be written in the following operator form 

(3.5) q a + + kA a u — 0, a — 1,2, 



G 
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du v-^ , 

(3-6) c — - J2 ^* a q a = /■ 

a— 1 

Thus, the system of equations governing the hyperbolic heat conduction does have 
the operator structure with conjugate operators. 

For hyperbolic heat conduction equation (|2.5I) the corresponding operator-differential 
equation has the form 

/q d 2 u du . r i df 



2 

(3.8) P=^2? Q , T> a u = A* a kA a , a =1,2. 

a=l 

Operator V, as well as its individual terms T> a , a = 1,2, is selfadjoint and positive 
definite in L 2 (il) on the set of functions satisfying boundary conditions (|2 . 3[) . 

We want to preserve the above operator structure of the differential model for the 
hyperbolic heat conduction after its approximation in space. For simplicity, we will 
consider the simplest difference approximations on uniform grids. In the considering 
problems it is natural to use for the scalar and vector unknowns staggered grids, where 
scalar variables and vector components employ their own grids. Such a technology is 
standard for problems of computational fluid dynamics |32j and electrodynamics |28] . 

The temperature is defined at the nodes of a uniform rectangular grid in Q: 

ui = {x I x = (xi, x 2 ), x a = i a h a , i a = 0, 1, N a , N a h a = l a , a = 1,2} 

and let w be a set of internal nodes (a) = uiUduj). The components of vector quantities 
are referred to the corresponding edges of the grid. We define 

Qi = {x I xi = (ii + 0.5)/ii, i\ = 0, 1, Ni — 1, x 2 = 12^2, «2 = 0, 1, N 2 } : 



u>2 = {x I xi = iihx, i\ = 0, 1, Nx, x 2 = («2 + 0.5)/i2, «2 = 0, 1, N 2 - 1} 

and Co a — u> a U du) a , a = 1,2. Component of heat flux q ai a — 1, 2 will be evaluated 
on the grid u a , a = 1, 2 (FigJBTT]). 

For grid functions y(x) = 0, x S 9w we define the Hilbert space if = £2(0;) with 
the scalar product and norm 

{y, w) = y(^)w(x)h 1 h 2 , \\y\\ = (y, y) 1/2 . 

Similarly, for the grid functions defined on grid u) a , a = 1,2, we define the Hilbert 
space H a , a = 1,2, where 

(y,w) a = Y y(x)w(x)h 1 h 2 , \\y\\ a = (y, yfj 2 , a = 1,2. 

Let us construct the grid analogs of differential operators A a , A* a , a = 1,2, 
defined above according to (|3.3[) . (|3.4I) . We will use the standard [18] central-difference 
approximations for derivatives in space. We set 

y(xi +0.5/ii,x 2 ) -y(xi -0.5/ii,x 2 ) 
(3.9) (Aij/)(x) = , xeui, 
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Fig. 3.1. 



uj(u), □ — wi(gi), O — o; 2 (g2). 



so that A\ : H —¥ Hi. Similarly, wc define A 2 : H ^ H 2 , where 
(3.10) 



y(a; 1 ,x 2 + 0.5/i 2 ) - J/(xi,x 2 - 0.5/i 2 ) 
(A 2 j/)(x) = , xew 2 . 



By the construction we have 

(3.11) A a tt = ^ Q u + 0(/i«), a = 1,2. 

Direct calculations verify that for the adjoint operators A* a : H a H, a = 1, 2 
we have the representation 

f«jio\ ^ n + 0.5/ii, x 2 ) - y(xi - 0.5hi,x 2 ) 

(3.12) (A x y)(x) = , xew, 



/ q1q n /„* w \ y(xi,x 2 + 0.5h 2 ) -y(xi,x 2 - 0.5/i 2 ) 

(3.13) (A 2 y) x = — i , xew. 

For sufficiently smooth functions u 

(3.14) A* a u = A* a u + 0{hl), a = 1,2 
holds. 

After approximation of system (|3.5|) . (|3.6|) in space we obtain the system of evo- 
lutionary equations 



(3.15) 



v— + kA a u — 0, x e w a , a = 1,2, 
at 
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(3.16) c ^-J2 A * a q a = f, xec. 

a— 1 

Similarly, equation ( ref (3.7)) is associated with the evolutionary equation 

d?u du Of 

(3.17) vc _ +c _ +Du = f + v _ xecj , 



(3.18) D = ^D a , D a =A* a kA a , a = 1,2. 

a=l 

Taking into account (|3~^1) . (|3~TU1) and (|3~T2"]l . (|3~T3f . for grid operators D a : H -> H, 
a = 1 , 2 we obtain 

(3.19) {Diy)(x) = — k{xi + 0.5hi, x 2 )- 



h\ \ " hi 



-k(xi — v.bhi,X2) ^ ) , x € cj. 



(3.20) (Z? 2 y)(x) = — \k{xx,x 2 + 0.5/i 2 ) 



/i2 V ^2 



t / n s2/(x) - y{x x ,x 2 - /i 2 )\ 
-fc(xi,x 2 - 0.5/i 2 ) : , xew. 

n-2 



Similarly to (l3~Til . (l3~14l) . we have [M[20] 

(3.21) D a u = V a u + 0{h 2 a ), a = 1,2 

in the class of sufficiently smooth coefficients k and functions it. In addition, in the 
space of grid functions H 

(3.22) D a = D* a , k S a E < D a < hA a E, 

x 4 . 2 irh a 4 2 nh a 

o a = — Sln "TTj 5 A Q = — cos , a = 1,2, 



2^ a 2/ Q 

where _E is the unit (identity) operator and fc(x) < ki, x G fi. 

4. Difference schemes for the hyperbolic heat conduction equation. We 

consider approximation in time for the approximate solution of differential-operator 
equation (|3.17[) . which is supplemented by the initial conditions 

du 

(4.1) m(x,0) = u (x), — (x,0) = vi(x), xew. 

at 

Let us define a uniform grid in time 

oj t = uj t U {T} = {t n = tit, n = 0,l,...,N, tN = T} 
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and denote y n = y(t n ), t n — nr. Standard three-level difference schemes with the 
second-order approximation in time will be considered. 

Equation (|3.7[) is approximated by the difference scheme with weights 

, s u n+1 - 2u n + u n - 1 
(4.2) vc - + c 



where 



2r 

+D{<ru n+1 + (1 - 2a)u n + au 11 - 1 ) = ip n , n = 1, 2, . . . , N - 1, 



ip n = f n + v J - 



2r 

with the corresponding initial data 

(4.3) w° = wo, = too. 



Scheme (|4.2[) . (|4.3[) which belongs to the class of three-level operator-difference 
schemes can be investigated on the basis of the Samarskii stability (correctness) theory 
of operator-difference schemes. Coincident necessary and sufficient conditions of the 
stability of these schemes in various norms are obtained in [T51 [H] . With this in mind, 
we give here only the simplest a priori estimates of stability with respect to the initial 
data and right hand side for scheme (I4.2p . (|4.3I) . 

Theorem 4.1. Difference scheme {4->ty i * s unconditionally stable at a > 

0.25 and for the finite- difference solution we have the estimate 

(4.4) S n+1 < S n + ^{c- 1 ^, <p n ), 



where 

(4.5) S n = [ [ vcE+ [ a- ~ ) t 2 D 



D- 



Proof. To prove this, we introduce the notation 



2 T 

Taking into account the identities 

u n = -(u n+1 + 2u n + u' 1 - 1 ) - -<u n+1 - 2u n + u"- 1 ), 
4 v 4 

au n+1 + (1 - 2a)u n + an"- 1 = u n + a(u n+1 - 2u n + u™" 1 ) 
we rewrite (14.21) in the form 



(4.6) [ vc E+[<j-~)t<D) 2 _ 1 + c ' - - +D - - - 
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We scalarly multiply in H this equation by 

2(C n+1 -C n )=r(Ti n+1 +Ti n )- 

This gives 

(4.7) (U+LJj^^H+pr 1 ,^ 1 ) 



vcE+[<T-±)r l D)ri n ,Ti n )-(DC,C) 



If g > 0.25 then value 

S n = ( (vcE + (a - T 2 D^j rt n ,rA + (DCC) 

defines the squared norm of the difference solution. With this notation we obtain the 
required estimate (14.51) . □ 

Estimate (|4. 5[) for the numerical solution is consistent with estimate (|2.12p for 
the solution of the differential problem. Using this estimate it is easy to prove in the 
standard enough way [T51 [05] that the difference solution converges to the exact one 
with truncation error 0(r 2 + h\ + (with the second order in time and space). 

5. Difference schemes for the hyperbolic heat conduction governed by 
the system of equations. For the approximate solution of the Cauchy problem for 
system fj3.5[) . (|3.6|) we use the simplest schemes with weights 

ri+l _ n 

(5.1) qZ {n) +v— — + kA a u a W = o, a = 1,2, 



(5.2) c ^ --£>«< B > = r +1/2 , n = ,l,...,J\T-l, 

a=l 

where a is a numerical parameter (weight), which is usually < a < 1. We have used 
the notation 

u a(n) = au n+l + ^ _ q a(n) = aq n+l + ^ _ ^ q n^ a=1}2 . 

For simplicity, we restrict ourselves to the same weight for all equations of system 
(j5~Tj) , (|iT2j) . Taking into account (f2T3l) we will supply (|5~Tj) . (|5"i>j) with the initial 
conditions 

(5.3) U ° = Uo , <Z°=5^ Q) , a = 1,2. 

We give the simplest estimates of stability for operator-difference scheme (|5.ip - (|5.3p . 
Estimate (|2.17p is used to guide us. 

Theorem 5.1. Difference scheme i5.1] )~ [5^3\) is unconditionally stable at a > 0.5 
and for the numerical solution the following estimate holds 

(5.4) G" +1 < exp (£\ G n +tT exp {c^P+^J^l 
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where 

2 

(5.5) G n = (cu n ,u n ) + uY,(k- 1 q^,q2)c 

a=l 

Proof. Scalarly multiply in H equation (|5.2[) by 2t , and each separate equa- 
tion (|5.1I) scalarly multiply in H a by 2rk~ l qa , a = 1,2 and sum them. Taking 
into account that 

2Tu a{n} = r{u n+1 + u n ) + (2(7 - 1)1 



T 

we obtain 

2 2 

(5.6) (cu n +\ u n+1 ) - (cu n , u n ) + vY, Q n a +1 ) a - v E( fc_1 9- C) a 



+ (2a-l)r 2 (c M — , — j+(2<7-l)r^X! ( fc r " g ° , Qa ~ 9a ) 



+2r^(fc- 1 ^("),^(")) a = 2r(/" +1 / 2 ,u^). 

a=l 

For terms in the right hand side of (|5.6|) we have 

2r(r +1 / 2 , u CT (")) = (2a - l)r 2 (V +1/2 , + T (f n+1 ' 2 , u n+1 + u n ). 

We restrict ourselves to schemes with a > 0.5 and use the estimates 

2 / m+1/2 U u_ 



(2a - l)r 2 / 



< (2a - l)r 



T(/ n+l/2 )U n+l +u „) < l. (c(/ +l +„»)) + I^( c -l r +l/2 ;r+ l/2 ); 



/ / n+1 i rc\ ^ ro+1 , n\\ ^ c\( n+1 n+l\ , o/ n n\ 
(c(m +11 ),(tl + U )) < 2(CU , U ) + 2(CU , U ). 

Substitution in (|5.6[) gives 

2 

(5-7) (l - 1) (c M " + V" +1 ) 
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< (l+I) (cu « lU »)+^(fc-i CgS ) a + !^ (i + ^r) (c-V" +1/2 ,r +1/2 ). 

a— 1 ^ ' 

Without loss of generality, we assume that 2r <T and therefore 

(l + ~) (l-^)"<exp(^). 

With this in mind, from (15.71) wc obtain timelevel-wise stability estimate (|5.4p . (|5 . 5[) . 
□ 

A priori estimate (|5.4j) is nothing but the grid analog of estimate (|2.17p and 
provides unconditional stability of the difference scheme with weights (|5.1I) . (|5.2[) 
under natural conditions a > 0.5. Considering the corresponding problem for the 
error [T5J [TH] , we prove the convergence of the solution of operator-difference problem 
(jB"3 ]) -(|0 ]) to the solution of differential-difference problem ([2~T]) . (|2~31) . (12TT31) at 
cr > 0.5 with order 0((2a — 1)t + r 2 ). If cr = 0.5, we have the second order of 
convergence with respect to r. 

The computational implementation of scheme (|5.1|) . (|5.2[) requires to solve the 
following grid problem at new time level n + 1 : 

(5.8) arql^ + vq n a +1 + a T kA a u n+1 = X n a , a = 1, 2, 



(5.9) cu n+1 -arJ2A* aq 2 +1 = r 

for given Xq, a = 1,2 and 0™. Substituting from equations (|5.8[) in equation 

(15.81). we obtain 



(5.10) (i/ + ar)cu n+1 + a 2 r 2 ^ ^A Q u n+1 = (i/ + ar)0" |ar^ A* xl- 

a— 1 a— 1 

Other components of the approximate solution are evaluated after solving grid prob- 
lem (|5.10|) via the explicit formulas of equations (|5.8[) . 

To preserve the second order of approximation, different grids in time are often 
employed for the individual components of the solution. The following scheme for 
system l|3.5p . (|3.6p provides an example 

n+l/2 n-1/2 n+1/2 n-1/2 

(5.11) Yv Y kA a u = 0, a = 1,2, 



(5.12) c " n+ ' T ^ -J2 A *^ +1/2 -f n+1/ ^ n = 0,l,...,N-l. 

a—l 

Such explicit schemes are widely used in computational practice. A detailed discus- 
sion of such schemes in application to problems of electrodynamics is presented, for 
example, in [55], with references to works of other authors. The main drawback of 
such schemes is connected with restrictions on the time step (conditional stability). 
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Considering equation ( ref (5.12)) at two time levels, we obtain the following 
equations 

2 n+1/2 . n-1/2 m+1/2 i fn-1/2 

E 4 . ga + ga _ J + J 

Q 2 ~ 2 ' 

Q = l 



n+1 o„Ti , „n-l 2 ™+l/2 o"- 1 / 2 fn+1/2 fn-1/2 

T 2 A " 7 ~~ 7 ■ 



U 

c 

a — 1 

Taking into account equation (|5.11[) . we derive 

. . u n+1 - 2u n + u"- 1 A 

(5.13) i/c ^ + c — + 2]i*M a « 11 

"~ a— 1 



jn+1/2 _|_ tn-1/2 jn+1/2 _ jn-1/2 
= '- h f . 

2 r 

Thus we have the explicit approximation of hyperbolic equation (|4.2|) with cr = 0. 
Unconditionally stable (at a > 0.25) scheme with weights 

+ V, U + y A * kAa{au n+l + 2a)u n + ^n-1) 

t z 2r ' 

a— 1 

jn+1/2 _|_ jn-1/2 jn+1/2 _ jn-1/2 

= 2 + " r 

is equivalent to the following scheme for system (|3.5|) , (|3.6[) , if in scheme (|5.1ip , (|5.12|l 
instead of (|5.11j) we use 

n+l/2 , n-1/2 n+1/2 n-1/2 

(5.14) + ^ 



+/b4 Q (cru" +1 + (1 - 2a)u n + aiS 1 - 1 ) = 0, a = 1, 2. 

Scheme (I5.12|) . f|5 . 14f) is not very convenient for the practical usage. Its main 
drawback results from the explicit coupling of equations for the temperature and heat 
fluxes. We must perform some preliminary work in order to obtain acceptable grid 
problems for evaluating the individual components of the solution at the new time 
level. 

Starting from explicit scheme (|5.11j) . (|5.12|) . we can construct unconditionally 
stable implicit schemes. We can do it in the most simple way using the Samarskii 
principle of regularization for difference schemes [181 119j , which is based on increasing 
the stability of a scheme via the perturbation of its operators. Stability of scheme 
(I5.11|) . (|5.12[) can be achieved in different ways. The most interesting possibility 
is connected with the multiplicative [22l [24] perturbation (increasing) of the time 
derivative operator or perturbation (decreasing) of the spatial variables operator for 
the individual equations of the system. 
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Consider the perturbation of equation (|5.12[) in detail. The implicit scheme can 
be written as 

(2 \ n+1 n ^ 

E + or 2 £ KA a cV= — ~ U - A *^ +1 ' 2 = f U+1/2 - 
a=l / a=l 

The perturbation has the order of C(r 2 ) and therefore regularized scheme (|5.1ip . 
(|5.15p remains in the class of schemes with the second order approximation. The 
operator at the time derivative is selfadjoint and positive definite. 

Theorem 5.2. Difference scheme i5.11\) , i5.15\) is unconditionally stable at 
vcqo > 0.25 and for the difference solution we have estimate j4-4^ where 

(5.16) S n = ((u<?-l*(E + a T 2 D)c^ 2 - T —B 



D u n +u n - 1 u n + u n ~ 1 



9 ' 



jn+l/2 _|_ yri-1/2 jn+l/2 _ jn-1/2 



Proof. From (|5 . 11 1) , (j5 . 15[) in the standard way we obtain the following three-level 
scheme 

1/9, 9 x 1/9 U U+1 -2u n + U Tl - 1 

(5.17) vc 1/2 (E + <jt 2 D)c 1/2 — 



+c 1/2 (E + <jt 2 D)c 1/2 - — ^ + Du n = y n , n = 1, 2, . . . , N - 1. 

2t 

Further investigation is conducted similarly to the proof of Theorem 14.11 Under the 
above restrictions on the weight we have 

vc 1 ^ 2 {E + ar 2 D)c 1 ^ 2 — — D > 
4 

and we associate the squared norm of the difference solution with S n . □ 

The numerical implementation of scheme (|5.1ip . (|5.15p is based on inversion of 
the same grid elliptic operator E + <tt 2 D, whereas the schemes with weights for the 
hyperbolic heat conduction equation (|4.2[) requires to invert c{y + t)E + gt 2 D. You 
can also obtain a grid analog of (|2.16[) . (|2.17p for difference scheme (|5.1ip . (|5.15|) . 
However, it seems difficult to proof same analog of Theorem 15. II in this. 

6. Splitting scheme for the hyperbolic heat conduction equation. The 

above considered unconditionally stable operator-difference schemes — (|4.2p for the 
hyperbolic heat equation and (|5.1ip , (|5. 15|) for the system of hyperbolic heat conduc- 
tion, respectively, — are not very convenient in the numerical implementation. We 
construct the additive schemes for problem (|2.5p , (|2.7p , where the transition to a new 
time level will be connected with the solution of more simple problems related to the 
inversion of individual operators A* a A ai a = 1,2, rather than their sum (operator 
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D in (I4.2p ). Taking into account the nature of operators A* a ,A ai a — 1,2, we are 
talking about locally one-dimensional schemes [T8] . 

We will focus on using regularized additive schemes of full approximation |23[ 125) . 
The principle of regularization of difference schemes is used traditionally widely |18j 
to construct stable difference schemes for the numerical solution of problems governed 
by partial differential equations. Due to small perturbations of the problem operators 
we can control the growth of the norm for the solution at the transition from one time 
level to another. 

The construction of unconditionally stable difference schemes via the principle of 
regularization is implemented as follows. For the initial problem there is constructed 
some simple difference scheme (producing difference scheme) which does not meet the 
necessary properties, ie the scheme is conditionally stable or even absolutely unstable. 
Then the quality of the difference scheme (its stability) is improved via perturbations 
of the difference scheme operators. 

It is natural to consider as the producing schemes the following explicit scheme 

u n+1 - 2u n +U 11 - 1 u^-u 11 - 1 

6.1 uc + c + Du n = <p r \ 

t z At 

which is complemented by initial conditions (|4.3p . The stability of this scheme (see 
(|4.3I) at a = 0) will be provided if the following inequality holds 

(6.2) R = vcE~- —D > 0. 

v ' 4 



In this case we have estimate (|4.4I) , in which 

an (u n -u n ~ x u n -u n - 1 \ / u n +u n - 1 u n +u n - 1 
(6.3) b = R , J + \D 



Taking into account (|3.22[) , from (|6.2j) we obtain the condition for stability of explicit 
scheme (|6.1|) 

~ fei(Ai + A 2 ) 1 2 

To increase the stability limit (increase operator i?), we can employ the perturbation 
of both the first term in R (ycE) and second one (D). 

In the case of perturbing the operator for the second time derivative we construct 
the regularized scheme by analogy with (|5.11[) . (|5.15|) : 

(6.4) vc^Qc 1 ' 2 - ~ 2 " + " + cV 2 Q C V 2 - ^— + Du n = <p n , 

t 2 2r 

where operator Q = Q* = E + 0(t 2 ). In the construction of additive schemes, we 
need to take into account the structure of the grid operator at new time level. Assume 
that 

(6.5) Q=(E+ ^t 2 A* 2 A 2 ) (E + <tt 2 A\Ai) (e + G -r 2 A\A 2 

so that Q > E + (tt 2 D. Direct calculations verify that for scheme (|6.4[) , (|6.5p instead 
of (16.21) we have 



(6.6) R = vc l ' 2 Qc 1 ' 2 - —D 
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and R > at vc a > 0.25. 

Theorem 6.1. Additive-difference scheme |6'.4| ) ; I 6. 5}) is unconditionally stable 
at vcqo > 0.25 and estimate j4-4\ ), \6.3\) , \6.b}) is valid for the difference solution. 

The second possibility of constructing unconditionally stable additive operator- 
difference schemes is connected with the perturbation of operator D in explicit scheme 
(|6.ip . Instead of operator D, which is defined according to (|3.18|) . we use 



2 
a=l 



«i.7l C=y]C a , C a = A* a (k- 1 E + ar 2 A a A* 

For these difference operators 



A n 



1,2. 



0"T Z 



1,2. 



Because of this, for 
(6.8) 

we have R > at i^coct > 0.5. 
Theorem 6.2. Additive 



R = vcE 



-C 



(6.9) 



2u n 



'erence schemes 



,n+l 



2r 



Cu n 



where operator C is defined according to \6. 7| ), is unconditionally stable at vcqCT > 0.5. 
Estimate J^.^| ) is frwe /or f/ie difference solution with 



(6.10) 5"= i?- — — + C- 



and i? corresponding to 116. 8\) . 

The main computational cost in the constructed splitting schemes results from 
the inversion of one-dimensional grid operators k~ 1 E + ot 1 A a A* a , a — 1,2. The 
potential advantage of additive scheme (|6.7|) . (|6.9|) in compare with scheme (|6.4p . 
(|6.5p is connected primarily with lower computational cost during the transition to a 
new time level. This advantage is more impressive for the three-dimensional problems 
(splitting in three directions). 

7. Additive schemes for the hyperbolic heat conduction governed by 

the system of equations. In the construction of splitting schemes for the system of 
equations governing the hyperbolic heat conduction, the theory and practice of using 
additive schemes for first order evolutionary equations will be employed. We treat 
the system of equations (|3.15[) , (|3.16|) as a single evolutionary equation for the vector 
u eee [tt 1 ,U2,Us] T = [qi,q 2 ,u] T : 

(7.1) B^+Au = f(t), 



where f = [0,0, f] T . For the elements of operator matrices B and A we have the 
representation 







" vk- 1 





" 




k^ 1 Ai 


(7.2) 


B = 





vk- 1 





A = 


k- 1 A 2 












c 




-A\ -A% 
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For the direct sum of spaces H = Hi © H 2 © H , we set 

p p 
(u,v) = ^2{u a ,v a ) a , ||u|| 2 = \\ U <*\\1- 

a—1 a—1 

In this case, A > in H and estimate (|2.16p . (|2.17[) can be rewritten as 

/ 

(7.3) (Bu(i), u(t)) < exp(t) (Bu(0), u(0)) + J exp(i - 6) (B _1 f (0), f (6»)) 2 d0. 

o 

To construct locally one-dimensional schemes for the Cauchy problem for (|7.ip . 
(|7.2[) . we use the additive representation of operator A in the form 

p 

(7.4) A=J2 A(a) - 

a=l 

The first variant of decomposition (|7.4[) corresponds to the selection p = 3 and 



(7.5) 



A« = 



Ax 


-A* 



A< 2 > = 




A 2 

-a; 



A (3) 



k- 1 
fc" 1 




Thus, we separate the individual terms with operators A a , A* a , a = 1,2. The main 
properties of these operators are connected with their non-negativity 

A( Q ) = -(A< Q ))*, a = 1, 2, A< 3 ) = (A< 3 ))* > 0. 

We can consider the splitting of (|7.4I) so that p — 2 where 



(7.6) A^ = 



0.5k- 1 





Ax 





0.5/e- 1 





-A\ 









A( 2 ) = 



o.sfc- 1 








0.5/fc- 1 
-AS 





A 2 




In this case A^ Q ) > 0, a = 1, 2 in space H. 

This non-negative property of operators A^ a \ a — 1,2, ...,p in splitting (|7.4j) 
allows to use for the approximate solution of the Cauchy problem for equation (|7.1| , 
(I7.2jl different classes of unconditionally stable additive operator-difference schemes 
[TS1|51]. With regard to our problem of the hyperbolic heat conduction, we can employ 
the schemes of second-order approximation in time. 

For the general case (in (17.41) p > 2) the standard additive schemes are based on 
the concept of summarized approximation. To construct the schemes of second order, 
we arrange computations via the algorithm 



1 



A (D _> I A (2) 



: A (p) I A bO ^ I A (p-U . . , _> I A (D 
! 2 2 2 
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The corresponding additive scheme of component-wise splitting seems like this: 

n+a/(2p) _ n+(a-l)/(2p) 1 _ 

(7.7) B- - h -A (a) (u" +a/(2p) + u "+ (a ~ 1)/(2p) ) 

r 2 

= S2 +1/ \ a=l,2,...,2p, 
where A< a ) = A^, a = 1, 2, A^ = A^ 1 ""), a = p + l,p + 2, 2p 

2p-2 

p+1/2 _ V"^ p+1/2 
a=l 

The proof of stability and convergence is conducted in the standard way, the 
technical details can be found, for example, in |241 . Scalarly multiplying the equations 
of scheme (|7.7p by y n + Q /( 2 P) + y"+(«- 1 )/( 2 p) ; we obtain the corresponding analog of 
a priori estimate (|7.3I) . 

Theorem 7.1. Additive operator- difference scheme of summarized approximation 
J7. 7| ) is unconditionally stable and approximates the system of equations |7. i| ) mi/i 
f/ie second order relative to r. 

The computational implementation of the considered additive schemes is much 
simpler than for schemes (|5.1I) . (|5.2p or (|5.1ip . f|5.15|) . To explain this fact, we consider, 
for example, the first step in scheme (17.71) with splitting (17.5[) . where 

(B + -A«) u" +1 / 6 = r" +1 / 2 

for given right hand side r n+1 / 2 . In the coordinate- wise form of this equation we have 
the system of equations 

n+l/6 . T . n+l/6 , n+1/2 n+l/6 , n+1/2 

J/Uj + —kAiy 3 = /c^ , = kr 2 , 

~."+V6 r A*i, n+1/6 - r n+1/2 
cu 3 ~ 2 11 ~~ 3 

Substituting u™ +1 ^ from the first equation into the third one, we obtain 

vcu\ l+1/& + ^AlkA lU n 3 +1/6 = vr n 3 +1/2 + ^-Alkr? +1/2 . 

Thus we must solve the one-dimensional grid problems with a single operator, which 
are connected with operators A\,A\. We have a similar realization for splitting (|7.6p . 

Among shortcomings of the constructed locally one-dimensional schemes (|7.7p 
it should be noted the lack of transparency (each individual equation does not ap- 
proximate the differential problem) as well as the relative difficulty of obtaining and 
studying schemes of increased approximation order. It is possible to construct for sys- 
tem (|3.15|) . (|3.16p another splitting schemes which belong to the class of regularized 
additive operator-difference schemes |24) . 

Regularized additive schemes can be constructed on the basis of scheme (|5.1ip . 
(|5T5j) . Instead of (|5T5)) we use the difference equation 

(7.8) c^Qc 1 ' 2 ^11^1 _ j2 a;c +1/2 - r +1/ '\ 
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where the factorized operator Q is defined according to (|6.5|) . Similar to Theorem l5.2[ 
we can prove the following statement involving estimate Q > E + ar 2 D. 

Theorem 7.2. Additive- difference scheme h5.11}) , \6.5\) . is unconditionally 

stable at vcqO > 0.25, and estimate |^.^[) holds for the difference solution, where 



S n =(( vc^Qc 1 ' 2 - -D 



D- 



Additive scheme (|5.11[) . (16.5[) . (|7.8|) is based on the perturbation of the operator at 
the time derivative in the last equation of system (|5.11l) . (|5.12[) (difference derivative 
of the temperature). It is interesting to consider the schemes with the perturbation 
of difference derivatives in time for heat fluxes. Instead of (|5.1ip we use the difference 
equations 

n+l/2 n— 1/2 

(7.9) {k ^E + ar 2 A a A* a )^ ±*i 



n+l/2 n— 1/2 

+is(k~ 1 + ar 2 A a A* a ) — ^ + A a u n — 0, a = 1,2. 

T 

Difference equations (|7.9[) can be written in the form 

n+l/2 , n-1/2 n+l/2 n-1/2 

q a + q a q a — Qa 

h V 



+ (k- 1 +<TT 2 A a A* a )~ 1 A a u n = 0, a = 1,2 

with treating them as the multiplicative perturbation of operators A a , a — 1, 2 (|5 . 1 1 [) . 

It is easy to see by direct calculations that additive scheme (|5.12l) (|7.9[) corresponds 
to scheme (jtTT)) ([Oil . 

Theorem 7.3. Additive-difference scheme \5.1'2\) , |7.ff[ ) is unconditionally stable 
at vcqO~ > 0.5, and estimate i6.10\) holds for the difference solution under the condition 
of setting operators C and R according to \6. 7| ) and W.8\) , respectively. 

In contrast to the scheme of summarized approximation (|7.7[) . operator-difference 
splitting scheme (|5.1ip . (|7.9j) has a clear and transparent structure. 
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